package ca.uol.aig.fftpack;
/**
  * @author Baoshe Zhang
  * @author Astronomical Instrument Group of University of Lethbridge.
*/
class RealDoubleFFT_Mixed
{

/*-------------------------------------------------
   radf2: Real FFT's forward processing of factor 2
  -------------------------------------------------*/
      void radf2(int ido, int l1, final double cc[], double ch[], 
                        final double wtable[], int offset)
      {
          int     i, k, ic;
          double  ti2, tr2;
          int iw1;
          iw1 = offset;

    	  for(k=0; k<l1; k++)
          {
	      ch[2*k*ido]=cc[k*ido]+cc[(k+l1)*ido];
	      ch[(2*k+1)*ido+ido-1]=cc[k*ido]-cc[(k+l1)*ido];
          }
          if(ido<2) return;
          if(ido !=2)
          {
	      for(k=0; k<l1; k++)
	      {
	          for(i=2; i<ido; i+=2)
	          {
		      ic=ido-i;
		      tr2 = wtable[i-2+iw1]*cc[i-1+(k+l1)*ido]
                           +wtable[i-1+iw1]*cc[i+(k+l1)*ido];
		      ti2 = wtable[i-2+iw1]*cc[i+(k+l1)*ido]
                           -wtable[i-1+iw1]*cc[i-1+(k+l1)*ido];
		      ch[i+2*k*ido]=cc[i+k*ido]+ti2;
		      ch[ic+(2*k+1)*ido]=ti2-cc[i+k*ido];
		      ch[i-1+2*k*ido]=cc[i-1+k*ido]+tr2;
		      ch[ic-1+(2*k+1)*ido]=cc[i-1+k*ido]-tr2;
	          }
	      }
	      if(ido%2==1)return;
          }
          for(k=0; k<l1; k++)
          {
	      ch[(2*k+1)*ido]=-cc[ido-1+(k+l1)*ido];
	      ch[ido-1+2*k*ido]=cc[ido-1+k*ido];
          }
      } 

/*-------------------------------------------------
   radb2: Real FFT's backward processing of factor 2
  -------------------------------------------------*/
      void radb2(int ido, int l1, final double cc[], double ch[], 
                        final double wtable[], int offset)
      {
          int     i, k, ic;
          double  ti2, tr2;
          int iw1 = offset;

          for(k=0; k<l1; k++)
          {
	      ch[k*ido]=cc[2*k*ido]+cc[ido-1+(2*k+1)*ido];
	      ch[(k+l1)*ido]=cc[2*k*ido]-cc[ido-1+(2*k+1)*ido];
          }
          if(ido<2) return;
          if(ido !=2)
          {
	      for(k=0; k<l1;++k)
	      {
	          for(i=2; i<ido; i+=2)
	          {
		      ic=ido-i;
		      ch[i-1+k*ido]=cc[i-1+2*k*ido]+cc[ic-1+(2*k+1)*ido];
		      tr2=cc[i-1+2*k*ido]-cc[ic-1+(2*k+1)*ido];
		      ch[i+k*ido]=cc[i+2*k*ido]-cc[ic+(2*k+1)*ido];
		      ti2=cc[i+(2*k)*ido]+cc[ic+(2*k+1)*ido];
		      ch[i-1+(k+l1)*ido]=wtable[i-2+iw1]*tr2-wtable[i-1+iw1]*ti2;
		      ch[i+(k+l1)*ido]=wtable[i-2+iw1]*ti2+wtable[i-1+iw1]*tr2;
	          }
	      }
	      if(ido%2==1) return;
          }
          for(k=0; k<l1; k++)
          {
	      ch[ido-1+k*ido]=2*cc[ido-1+2*k*ido];
	      ch[ido-1+(k+l1)*ido]=-2*cc[(2*k+1)*ido];
          }
     }

/*-------------------------------------------------
   radf3: Real FFT's forward processing of factor 3 
  -------------------------------------------------*/
     void radf3(int ido, int l1, final double cc[], double ch[], 
                       final double wtable[], int offset)
     {
          final double taur=-0.5D;
          final double taui=0.866025403784439D;
          int     i, k, ic;
          double  ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
          int iw1, iw2;
          iw1 = offset;
          iw2 = iw1 + ido;

          for(k=0; k<l1; k++)
          {
	      cr2=cc[(k+l1)*ido]+cc[(k+2*l1)*ido];
	      ch[3*k*ido]=cc[k*ido]+cr2;
	      ch[(3*k+2)*ido]=taui*(cc[(k+l1*2)*ido]-cc[(k+l1)*ido]);
	      ch[ido-1+(3*k+1)*ido]=cc[k*ido]+taur*cr2;
          }
          if(ido==1) return;
          for(k=0; k<l1; k++)
          {
	      for(i=2; i<ido; i+=2)
	      {
	          ic=ido-i;
	          dr2 = wtable[i-2+iw1]*cc[i-1+(k+l1)*ido]
                       +wtable[i-1+iw1]*cc[i+(k+l1)*ido];
	          di2 = wtable[i-2+iw1]*cc[i+(k+l1)*ido]
                       -wtable[i-1+iw1]*cc[i-1+(k+l1)*ido];
	          dr3 = wtable[i-2+iw2]*cc[i-1+(k+l1*2)*ido]
                       +wtable[i-1+iw2]*cc[i+(k+l1*2)*ido];
	          di3 = wtable[i-2+iw2]*cc[i+(k+l1*2)*ido]
                       -wtable[i-1+iw2]*cc[i-1+(k+l1*2)*ido];
	          cr2 = dr2+dr3;
	          ci2 = di2+di3;
	          ch[i-1+3*k*ido]=cc[i-1+k*ido]+cr2;
	          ch[i+3*k*ido]=cc[i+k*ido]+ci2;
	          tr2=cc[i-1+k*ido]+taur*cr2;
	          ti2=cc[i+k*ido]+taur*ci2;
	          tr3=taui*(di2-di3);
	          ti3=taui*(dr3-dr2);
	          ch[i-1+(3*k+2)*ido]=tr2+tr3;
	          ch[ic-1+(3*k+1)*ido]=tr2-tr3;
	          ch[i+(3*k+2)*ido]=ti2+ti3;
	          ch[ic+(3*k+1)*ido]=ti3-ti2;
	      }
          }
     } 

/*-------------------------------------------------
   radb3: Real FFT's backward processing of factor 3
  -------------------------------------------------*/
     void radb3(int ido, int l1, final double cc[], double ch[], 
                       final double wtable[], int offset)
     {
          final double taur=-0.5D;
          final double taui=0.866025403784439D;
          int i, k, ic;
          double  ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
          int iw1, iw2;
          iw1 = offset;
          iw2 = iw1 + ido;

          for(k=0; k<l1; k++)
          {
	      tr2=2*cc[ido-1+(3*k+1)*ido];
	      cr2=cc[3*k*ido]+taur*tr2;
	      ch[k*ido]=cc[3*k*ido]+tr2;
	      ci3=2*taui*cc[(3*k+2)*ido];
	      ch[(k+l1)*ido]=cr2-ci3;
	      ch[(k+2*l1)*ido]=cr2+ci3;
          }
          if(ido==1) return;
          for(k=0; k<l1; k++)
          {
	      for(i=2; i<ido; i+=2)
	      {
	          ic=ido-i;
	          tr2=cc[i-1+(3*k+2)*ido]+cc[ic-1+(3*k+1)*ido];
	          cr2=cc[i-1+3*k*ido]+taur*tr2;
	          ch[i-1+k*ido]=cc[i-1+3*k*ido]+tr2;
	          ti2=cc[i+(3*k+2)*ido]-cc[ic+(3*k+1)*ido];
	          ci2=cc[i+3*k*ido]+taur*ti2;
	          ch[i+k*ido]=cc[i+3*k*ido]+ti2;
	          cr3=taui*(cc[i-1+(3*k+2)*ido]-cc[ic-1+(3*k+1)*ido]);
	          ci3=taui*(cc[i+(3*k+2)*ido]+cc[ic+(3*k+1)*ido]);
	          dr2=cr2-ci3;
	          dr3=cr2+ci3;
	          di2=ci2+cr3;
	          di3=ci2-cr3;
	          ch[i-1+(k+l1)*ido] = wtable[i-2+iw1]*dr2
                                      -wtable[i-1+iw1]*di2;
	          ch[i+(k+l1)*ido] = wtable[i-2+iw1]*di2
                                    +wtable[i-1+iw1]*dr2;
	          ch[i-1+(k+2*l1)*ido] = wtable[i-2+iw2]*dr3
                                        -wtable[i-1+iw2]*di3;
	          ch[i+(k+2*l1)*ido] = wtable[i-2+iw2]*di3
                                      +wtable[i-1+iw2]*dr3;
	      }
          }
     } 

/*-------------------------------------------------
   radf4: Real FFT's forward processing of factor 4
  -------------------------------------------------*/
     void radf4(int ido, int l1, final double cc[], double ch[], 
                       final double wtable[], int offset)
     {
          final double hsqt2=0.7071067811865475D;
          int i, k, ic;
          double  ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
          int iw1, iw2, iw3;
          iw1 = offset;
          iw2 = offset + ido;
          iw3 = iw2 + ido;
          for(k=0; k<l1; k++)
          {
	      tr1=cc[(k+l1)*ido]+cc[(k+3*l1)*ido];
	      tr2=cc[k*ido]+cc[(k+2*l1)*ido];
	      ch[4*k*ido]=tr1+tr2;
	      ch[ido-1+(4*k+3)*ido]=tr2-tr1;
	      ch[ido-1+(4*k+1)*ido]=cc[k*ido]-cc[(k+2*l1)*ido];
	      ch[(4*k+2)*ido]=cc[(k+3*l1)*ido]-cc[(k+l1)*ido];
          }
          if(ido<2) return;
          if(ido !=2)
          {
	      for(k=0; k<l1; k++)
	      {
	          for(i=2; i<ido; i+=2)
	          {
			ic=ido-i;
			cr2 = wtable[i-2+iw1]*cc[i-1+(k+l1)*ido]
                             +wtable[i-1+iw1]*cc[i+(k+l1)*ido];
			ci2 = wtable[i-2+iw1]*cc[i+(k+l1)*ido]
                             -wtable[i-1+iw1]*cc[i-1+(k+l1)*ido];
			cr3 = wtable[i-2+iw2]*cc[i-1+(k+2*l1)*ido]
                             +wtable[i-1+iw2]*cc[i+(k+2*l1)*ido];
			ci3 = wtable[i-2+iw2]*cc[i+(k+2*l1)*ido]
                             -wtable[i-1+iw2]*cc[i-1+(k+2*l1)*ido];
			cr4 = wtable[i-2+iw3]*cc[i-1+(k+3*l1)*ido]
                             +wtable[i-1+iw3]*cc[i+(k+3*l1)*ido];
			ci4 = wtable[i-2+iw3]*cc[i+(k+3*l1)*ido]
                             -wtable[i-1+iw3]*cc[i-1+(k+3*l1)*ido];
			tr1=cr2+cr4;
			tr4=cr4-cr2;
			ti1=ci2+ci4;
			ti4=ci2-ci4;
			ti2=cc[i+k*ido]+ci3;
			ti3=cc[i+k*ido]-ci3;
			tr2=cc[i-1+k*ido]+cr3;
			tr3=cc[i-1+k*ido]-cr3;
			ch[i-1+4*k*ido]=tr1+tr2;
			ch[ic-1+(4*k+3)*ido]=tr2-tr1;
			ch[i+4*k*ido]=ti1+ti2;
			ch[ic+(4*k+3)*ido]=ti1-ti2;
			ch[i-1+(4*k+2)*ido]=ti4+tr3;
			ch[ic-1+(4*k+1)*ido]=tr3-ti4;
			ch[i+(4*k+2)*ido]=tr4+ti3;
			ch[ic+(4*k+1)*ido]=tr4-ti3;
	          }
	      }
	      if(ido%2==1) return;
          }
          for(k=0; k<l1; k++)
          {
	      ti1=-hsqt2*(cc[ido-1+(k+l1)*ido]+cc[ido-1+(k+3*l1)*ido]);
              tr1=hsqt2*(cc[ido-1+(k+l1)*ido]-cc[ido-1+(k+3*l1)*ido]);
              ch[ido-1+4*k*ido]=tr1+cc[ido-1+k*ido];
	      ch[ido-1+(4*k+2)*ido]=cc[ido-1+k*ido]-tr1;
	      ch[(4*k+1)*ido]=ti1-cc[ido-1+(k+2*l1)*ido];
	      ch[(4*k+3)*ido]=ti1+cc[ido-1+(k+2*l1)*ido];
          }
     } 

/*-------------------------------------------------
   radb4: Real FFT's backward processing of factor 4
  -------------------------------------------------*/
     void radb4(int ido, int l1, final double cc[], double ch[], 
                       final double wtable[], int offset)
     {
         final double sqrt2=1.414213562373095D;
         int i, k, ic;
         double  ci2, ci3, ci4, cr2, cr3, cr4; 
         double  ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
         int iw1, iw2, iw3;
         iw1 = offset;
         iw2 = iw1 + ido;
         iw3 = iw2 + ido;
         
         for(k=0; k<l1; k++)
         {
	     tr1=cc[4*k*ido]-cc[ido-1+(4*k+3)*ido];
	     tr2=cc[4*k*ido]+cc[ido-1+(4*k+3)*ido];
	     tr3=cc[ido-1+(4*k+1)*ido]+cc[ido-1+(4*k+1)*ido];
	     tr4=cc[(4*k+2)*ido]+cc[(4*k+2)*ido];
	     ch[k*ido]=tr2+tr3;
	     ch[(k+l1)*ido]=tr1-tr4;
	     ch[(k+2*l1)*ido]=tr2-tr3;
	     ch[(k+3*l1)*ido]=tr1+tr4;
         }
         if(ido<2) return;
         if(ido !=2)
         {
	     for(k=0; k<l1;++k)
	     {
	         for(i=2; i<ido; i+=2)
	         {
			ic=ido-i;
			ti1=cc[i+4*k*ido]+cc[ic+(4*k+3)*ido];
			ti2=cc[i+4*k*ido]-cc[ic+(4*k+3)*ido];
			ti3=cc[i+(4*k+2)*ido]-cc[ic+(4*k+1)*ido];
			tr4=cc[i+(4*k+2)*ido]+cc[ic+(4*k+1)*ido];
			tr1=cc[i-1+4*k*ido]-cc[ic-1+(4*k+3)*ido];
			tr2=cc[i-1+4*k*ido]+cc[ic-1+(4*k+3)*ido];
			ti4=cc[i-1+(4*k+2)*ido]-cc[ic-1+(4*k+1)*ido];
			tr3=cc[i-1+(4*k+2)*ido]+cc[ic-1+(4*k+1)*ido];
			ch[i-1+k*ido]=tr2+tr3;
			cr3=tr2-tr3;
			ch[i+k*ido]=ti2+ti3;
			ci3=ti2-ti3;
			cr2=tr1-tr4;
			cr4=tr1+tr4;
			ci2=ti1+ti4;
			ci4=ti1-ti4;
			ch[i-1+(k+l1)*ido] = wtable[i-2+iw1]*cr2
                                            -wtable[i-1+iw1]*ci2;
			ch[i+(k+l1)*ido] = wtable[i-2+iw1]*ci2
                                          +wtable[i-1+iw1]*cr2;
			ch[i-1+(k+2*l1)*ido] = wtable[i-2+iw2]*cr3
                                              -wtable[i-1+iw2]*ci3;
			ch[i+(k+2*l1)*ido] = wtable[i-2+iw2]*ci3
                                            +wtable[i-1+iw2]*cr3;
			ch[i-1+(k+3*l1)*ido] = wtable[i-2+iw3]*cr4
                                              -wtable[i-1+iw3]*ci4;
			ch[i+(k+3*l1)*ido] = wtable[i-2+iw3]*ci4
                                            +wtable[i-1+iw3]*cr4;
	         }
	     }
	     if(ido%2==1) return;
         }
         for(k=0; k<l1; k++)
         {
	     ti1=cc[(4*k+1)*ido]+cc[(4*k+3)*ido];
	     ti2=cc[(4*k+3)*ido]-cc[(4*k+1)*ido];
	     tr1=cc[ido-1+4*k*ido]-cc[ido-1+(4*k+2)*ido];
	     tr2=cc[ido-1+4*k*ido]+cc[ido-1+(4*k+2)*ido];
	     ch[ido-1+k*ido]=tr2+tr2;
	     ch[ido-1+(k+l1)*ido]=sqrt2*(tr1-ti1);
	     ch[ido-1+(k+2*l1)*ido]=ti2+ti2;
	     ch[ido-1+(k+3*l1)*ido]=-sqrt2*(tr1+ti1);
         }
     } 

/*-------------------------------------------------
   radf5: Real FFT's forward processing of factor 5
  -------------------------------------------------*/
     void radf5(int ido, int l1, final double cc[], double ch[], 
                       final double wtable[], int offset)
     {
         final double tr11=0.309016994374947D;
         final double ti11=0.951056516295154D;
         final double tr12=-0.809016994374947D;
         final double ti12=0.587785252292473D;
         int     i, k, ic;
         double  ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3,
                 dr4, dr5, cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
         int iw1, iw2, iw3, iw4;
         iw1 = offset;
         iw2 = iw1 + ido;
         iw3 = iw2 + ido;
         iw4 = iw3 + ido;

         for(k=0; k<l1; k++)
         {
	     cr2=cc[(k+4*l1)*ido]+cc[(k+l1)*ido];
	     ci5=cc[(k+4*l1)*ido]-cc[(k+l1)*ido];
	     cr3=cc[(k+3*l1)*ido]+cc[(k+2*l1)*ido];
	     ci4=cc[(k+3*l1)*ido]-cc[(k+2*l1)*ido];
	     ch[5*k*ido]=cc[k*ido]+cr2+cr3;
	     ch[ido-1+(5*k+1)*ido]=cc[k*ido]+tr11*cr2+tr12*cr3;
	     ch[(5*k+2)*ido]=ti11*ci5+ti12*ci4;
	     ch[ido-1+(5*k+3)*ido]=cc[k*ido]+tr12*cr2+tr11*cr3;
	     ch[(5*k+4)*ido]=ti12*ci5-ti11*ci4;
         }
         if(ido==1) return;
         for(k=0; k<l1;++k)
         {
	     for(i=2; i<ido; i+=2)
	     {
	            ic=ido-i;
		    dr2 = wtable[i-2+iw1]*cc[i-1+(k+l1)*ido]
                         +wtable[i-1+iw1]*cc[i+(k+l1)*ido];
		    di2 = wtable[i-2+iw1]*cc[i+(k+l1)*ido]
                         -wtable[i-1+iw1]*cc[i-1+(k+l1)*ido];
		    dr3 = wtable[i-2+iw2]*cc[i-1+(k+2*l1)*ido]
                         +wtable[i-1+iw2]*cc[i+(k+2*l1)*ido];
		    di3 = wtable[i-2+iw2]*cc[i+(k+2*l1)*ido]
                         -wtable[i-1+iw2]*cc[i-1+(k+2*l1)*ido];
		    dr4 = wtable[i-2+iw3]*cc[i-1+(k+3*l1)*ido]
                         +wtable[i-1+iw3]*cc[i+(k+3*l1)*ido];
		    di4 = wtable[i-2+iw3]*cc[i+(k+3*l1)*ido]
                         -wtable[i-1+iw3]*cc[i-1+(k+3*l1)*ido];
		    dr5 = wtable[i-2+iw4]*cc[i-1+(k+4*l1)*ido]
                         +wtable[i-1+iw4]*cc[i+(k+4*l1)*ido];
		    di5 = wtable[i-2+iw4]*cc[i+(k+4*l1)*ido]
                         -wtable[i-1+iw4]*cc[i-1+(k+4*l1)*ido];
		    cr2=dr2+dr5;
		    ci5=dr5-dr2;
		    cr5=di2-di5;
		    ci2=di2+di5;
		    cr3=dr3+dr4;
		    ci4=dr4-dr3;
		    cr4=di3-di4;
		    ci3=di3+di4;
		    ch[i-1+5*k*ido]=cc[i-1+k*ido]+cr2+cr3;
		    ch[i+5*k*ido]=cc[i+k*ido]+ci2+ci3;
		    tr2=cc[i-1+k*ido]+tr11*cr2+tr12*cr3;
		    ti2=cc[i+k*ido]+tr11*ci2+tr12*ci3;
		    tr3=cc[i-1+k*ido]+tr12*cr2+tr11*cr3;
		    ti3=cc[i+k*ido]+tr12*ci2+tr11*ci3;
		    tr5=ti11*cr5+ti12*cr4;
		    ti5=ti11*ci5+ti12*ci4;
		    tr4=ti12*cr5-ti11*cr4;
		    ti4=ti12*ci5-ti11*ci4;
		    ch[i-1+(5*k+2)*ido]=tr2+tr5;
		    ch[ic-1+(5*k+1)*ido]=tr2-tr5;
		    ch[i+(5*k+2)*ido]=ti2+ti5;
		    ch[ic+(5*k+1)*ido]=ti5-ti2;
		    ch[i-1+(5*k+4)*ido]=tr3+tr4;
		    ch[ic-1+(5*k+3)*ido]=tr3-tr4;
		    ch[i+(5*k+4)*ido]=ti3+ti4;
		    ch[ic+(5*k+3)*ido]=ti4-ti3;
	     }
         }
     } 

/*-------------------------------------------------
   radb5: Real FFT's backward processing of factor 5
  -------------------------------------------------*/
     void radb5(int ido, int l1, final double cc[], double ch[], 
                       final double wtable[], int offset)
     {
         final double tr11=0.309016994374947D;
         final double ti11=0.951056516295154D;
         final double tr12=-0.809016994374947D;
         final double ti12=0.587785252292473D;
         int     i, k, ic;
         double  ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4,
                 ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
         int iw1, iw2, iw3, iw4;
         iw1 = offset;
         iw2 = iw1 + ido;
         iw3 = iw2 + ido;
         iw4 = iw3 + ido;

         for(k=0; k<l1; k++)
         {
		ti5=2*cc[(5*k+2)*ido];
		ti4=2*cc[(5*k+4)*ido];
		tr2=2*cc[ido-1+(5*k+1)*ido];
		tr3=2*cc[ido-1+(5*k+3)*ido];
		ch[k*ido]=cc[5*k*ido]+tr2+tr3;
		cr2=cc[5*k*ido]+tr11*tr2+tr12*tr3;
		cr3=cc[5*k*ido]+tr12*tr2+tr11*tr3;
		ci5=ti11*ti5+ti12*ti4;
		ci4=ti12*ti5-ti11*ti4;
		ch[(k+l1)*ido]=cr2-ci5;
		ch[(k+2*l1)*ido]=cr3-ci4;
		ch[(k+3*l1)*ido]=cr3+ci4;
		ch[(k+4*l1)*ido]=cr2+ci5;
         }
         if(ido==1) return;
         for(k=0; k<l1;++k)
         {
	     for(i=2; i<ido; i+=2)
	     {
		    ic=ido-i;
		    ti5=cc[i+(5*k+2)*ido]+cc[ic+(5*k+1)*ido];
		    ti2=cc[i+(5*k+2)*ido]-cc[ic+(5*k+1)*ido];
		    ti4=cc[i+(5*k+4)*ido]+cc[ic+(5*k+3)*ido];
		    ti3=cc[i+(5*k+4)*ido]-cc[ic+(5*k+3)*ido];
		    tr5=cc[i-1+(5*k+2)*ido]-cc[ic-1+(5*k+1)*ido];
		    tr2=cc[i-1+(5*k+2)*ido]+cc[ic-1+(5*k+1)*ido];
		    tr4=cc[i-1+(5*k+4)*ido]-cc[ic-1+(5*k+3)*ido];
		    tr3=cc[i-1+(5*k+4)*ido]+cc[ic-1+(5*k+3)*ido];
		    ch[i-1+k*ido]=cc[i-1+5*k*ido]+tr2+tr3;
		    ch[i+k*ido]=cc[i+5*k*ido]+ti2+ti3;
		    cr2=cc[i-1+5*k*ido]+tr11*tr2+tr12*tr3;

		    ci2=cc[i+5*k*ido]+tr11*ti2+tr12*ti3;
		    cr3=cc[i-1+5*k*ido]+tr12*tr2+tr11*tr3;

		    ci3=cc[i+5*k*ido]+tr12*ti2+tr11*ti3;
		    cr5=ti11*tr5+ti12*tr4;
		    ci5=ti11*ti5+ti12*ti4;
		    cr4=ti12*tr5-ti11*tr4;
		    ci4=ti12*ti5-ti11*ti4;
		    dr3=cr3-ci4;
		    dr4=cr3+ci4;
		    di3=ci3+cr4;
		    di4=ci3-cr4;
		    dr5=cr2+ci5;
		    dr2=cr2-ci5;
		    di5=ci2-cr5;
		    di2=ci2+cr5;
		    ch[i-1+(k+l1)*ido] = wtable[i-2+iw1]*dr2
                                        -wtable[i-1+iw1]*di2;
		    ch[i+(k+l1)*ido] = wtable[i-2+iw1]*di2
                                      +wtable[i-1+iw1]*dr2;
		    ch[i-1+(k+2*l1)*ido] = wtable[i-2+iw2]*dr3
                                          -wtable[i-1+iw2]*di3;
		    ch[i+(k+2*l1)*ido] = wtable[i-2+iw2]*di3
                                        +wtable[i-1+iw2]*dr3;
		    ch[i-1+(k+3*l1)*ido] = wtable[i-2+iw3]*dr4
                                          -wtable[i-1+iw3]*di4;
		    ch[i+(k+3*l1)*ido] = wtable[i-2+iw3]*di4
                                        +wtable[i-1+iw3]*dr4;
		    ch[i-1+(k+4*l1)*ido] = wtable[i-2+iw4]*dr5
                                          -wtable[i-1+iw4]*di5;
		    ch[i+(k+4*l1)*ido] = wtable[i-2+iw4]*di5
                                        +wtable[i-1+iw4]*dr5;
	     }
        }
     } 

/*---------------------------------------------------------
   radfg: Real FFT's forward processing of general factor
  --------------------------------------------------------*/
     void radfg(int ido, int ip, int l1, int idl1, double cc[], 
                       double c1[], double c2[], double ch[], double ch2[], 
                       final double wtable[], int offset)
     {
          final double twopi=2.0D*Math.PI; //6.28318530717959;
          int     idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;
          double  dc2, ai1, ai2, ar1, ar2, ds2, dcp, arg, dsp, ar1h, ar2h;
          int iw1 = offset;

          arg=twopi / (double)ip;
          dcp=Math.cos(arg);
          dsp=Math.sin(arg);
          ipph=(ip+1)/ 2;
          nbd=(ido-1)/ 2;
          if(ido !=1)
          {
	      for(ik=0; ik<idl1; ik++) ch2[ik]=c2[ik];
	      for(j=1; j<ip; j++)
	         for(k=0; k<l1; k++)
		     ch[(k+j*l1)*ido]=c1[(k+j*l1)*ido];
	      if(nbd<=l1)
	      {
	          is=-ido;
	          for(j=1; j<ip; j++)
	          {
		      is+=ido;
		      idij=is-1;
		      for(i=2; i<ido; i+=2)
		      {
		          idij+=2;
		          for(k=0; k<l1; k++)
		          {
			       ch[i-1+(k+j*l1)*ido]=
			            wtable[idij-1+iw1]*c1[i-1+(k+j*l1)*ido]
                                   +wtable[idij+iw1]*c1[i+(k+j*l1)*ido];
			       ch[i+(k+j*l1)*ido]=
			            wtable[idij-1+iw1]*c1[i+(k+j*l1)*ido]
                                   -wtable[idij+iw1]*c1[i-1+(k+j*l1)*ido];
		          }
		      }
	          }
	      }
	      else
	      {
	          is=-ido;
	          for(j=1; j<ip; j++)
	          {
		      is+=ido;
		      for(k=0; k<l1; k++)
		      {
		          idij=is-1;
		          for(i=2; i<ido; i+=2)
		          {
			      idij+=2;
			      ch[i-1+(k+j*l1)*ido]=
			           wtable[idij-1+iw1]*c1[i-1+(k+j*l1)*ido]
                                  +wtable[idij+iw1]*c1[i+(k+j*l1)*ido];
			      ch[i+(k+j*l1)*ido]=
			           wtable[idij-1+iw1]*c1[i+(k+j*l1)*ido]
                                  -wtable[idij+iw1]*c1[i-1+(k+j*l1)*ido];
		          }
		      }
	          }
	      }
	      if(nbd>=l1)
	      {
	          for(j=1; j<ipph; j++)
	          {
		      jc=ip-j;
		      for(k=0; k<l1; k++)
		      {
		          for(i=2; i<ido; i+=2)
		          {
			      c1[i-1+(k+j*l1)*ido]=ch[i-1+(k+j*l1)*ido]+ch[i-1+(k+jc*l1)*ido];
			      c1[i-1+(k+jc*l1)*ido]=ch[i+(k+j*l1)*ido]-ch[i+(k+jc*l1)*ido];
			      c1[i+(k+j*l1)*ido]=ch[i+(k+j*l1)*ido]+ch[i+(k+jc*l1)*ido];
			      c1[i+(k+jc*l1)*ido]=ch[i-1+(k+jc*l1)*ido]-ch[i-1+(k+j*l1)*ido];
		          }
		      }
	          }
	      }
	      else
	      {
	          for(j=1; j<ipph; j++)
	          {
		      jc=ip-j;
		      for(i=2; i<ido; i+=2)
		      {
		          for(k=0; k<l1; k++)
		          {
			      c1[i-1+(k+j*l1)*ido]=
			          ch[i-1+(k+j*l1)*ido]+ch[i-1+(k+jc*l1)*ido];
			      c1[i-1+(k+jc*l1)*ido]=ch[i+(k+j*l1)*ido]-ch[i+(k+jc*l1)*ido];
			      c1[i+(k+j*l1)*ido]=ch[i+(k+j*l1)*ido]+ch[i+(k+jc*l1)*ido];
			      c1[i+(k+jc*l1)*ido]=ch[i-1+(k+jc*l1)*ido]-ch[i-1+(k+j*l1)*ido];
		          }
		      }
	          }
	      }
          }
          else
          {				
	      for(ik=0; ik<idl1; ik++) c2[ik]=ch2[ik];
          }
          for(j=1; j<ipph; j++)
          {
	      jc=ip-j;
	      for(k=0; k<l1; k++)
	      {
	          c1[(k+j*l1)*ido]=ch[(k+j*l1)*ido]+ch[(k+jc*l1)*ido];
	          c1[(k+jc*l1)*ido]=ch[(k+jc*l1)*ido]-ch[(k+j*l1)*ido];
	      }
          }

          ar1=1;
          ai1=0;
          for(l=1; l<ipph; l++)
          {
	      lc=ip-l;
	      ar1h=dcp*ar1-dsp*ai1;
	      ai1=dcp*ai1+dsp*ar1;
	      ar1=ar1h;
	      for(ik=0; ik<idl1; ik++)
	      {
	          ch2[ik+l*idl1]=c2[ik]+ar1*c2[ik+idl1];
	          ch2[ik+lc*idl1]=ai1*c2[ik+(ip-1)*idl1];
	      }
	      dc2=ar1;
	      ds2=ai1;
	      ar2=ar1;
	      ai2=ai1;
	      for(j=2; j<ipph; j++)
	      {
	          jc=ip-j;
	          ar2h=dc2*ar2-ds2*ai2;
	          ai2=dc2*ai2+ds2*ar2;
	          ar2=ar2h;
	          for(ik=0; ik<idl1; ik++)
	          {
		      ch2[ik+l*idl1]+=ar2*c2[ik+j*idl1];
		      ch2[ik+lc*idl1]+=ai2*c2[ik+jc*idl1];
	          }
	      }
          }
          for(j=1; j<ipph; j++)
	     for(ik=0; ik<idl1; ik++)
	        ch2[ik]+=c2[ik+j*idl1];

          if(ido>=l1)
          {
	     for(k=0; k<l1; k++)
	     {
	        for(i=0; i<ido; i++)
	        {
		    cc[i+k*ip*ido]=ch[i+k*ido];
	        }
	     }
          }
          else
          {
	      for(i=0; i<ido; i++)
	      {
	          for(k=0; k<l1; k++)
	          {
		      cc[i+k*ip*ido]=ch[i+k*ido];
	          }
	      }
          }
          for(j=1; j<ipph; j++)
          {
	      jc=ip-j;
	      j2=2*j;
	      for(k=0; k<l1; k++)
	      {
	          cc[ido-1+(j2-1+k*ip)*ido]=ch[(k+j*l1)*ido];
	          cc[(j2+k*ip)*ido]=ch[(k+jc*l1)*ido];
	      }
          }
          if(ido==1) return;
          if(nbd>=l1)
          {
	      for(j=1; j<ipph; j++)
	      {
	          jc=ip-j;
	          j2=2*j;
	          for(k=0; k<l1; k++)
	          {
		      for(i=2; i<ido; i+=2)
		      {
		          ic=ido-i;
		          cc[i-1+(j2+k*ip)*ido]=ch[i-1+(k+j*l1)*ido]+ch[i-1+(k+jc*l1)*ido];
		          cc[ic-1+(j2-1+k*ip)*ido]=ch[i-1+(k+j*l1)*ido]-ch[i-1+(k+jc*l1)*ido];
		          cc[i+(j2+k*ip)*ido]=ch[i+(k+j*l1)*ido]+ch[i+(k+jc*l1)*ido];
		          cc[ic+(j2-1+k*ip)*ido]=ch[i+(k+jc*l1)*ido]-ch[i+(k+j*l1)*ido];
		      }
	          }
	      }
          }
          else
          {
	      for(j=1; j<ipph; j++)
	      {
	          jc=ip-j;
	          j2=2*j;
	          for(i=2; i<ido; i+=2)
	          {
		      ic=ido-i;
		      for(k=0; k<l1; k++)
		      {
		          cc[i-1+(j2+k*ip)*ido]=ch[i-1+(k+j*l1)*ido]+ch[i-1+(k+jc*l1)*ido];
		          cc[ic-1+(j2-1+k*ip)*ido]=ch[i-1+(k+j*l1)*ido]-ch[i-1+(k+jc*l1)*ido];
		          cc[i+(j2+k*ip)*ido]=ch[i+(k+j*l1)*ido]+ch[i+(k+jc*l1)*ido];
		          cc[ic+(j2-1+k*ip)*ido]=ch[i+(k+jc*l1)*ido]-ch[i+(k+j*l1)*ido];
		      }
	          }
	      }
          }
     } 

/*---------------------------------------------------------
   radbg: Real FFT's backward processing of general factor
  --------------------------------------------------------*/
     void radbg(int ido, int ip, int l1, int idl1, double cc[], double c1[], 
                       double c2[], double ch[], double ch2[], final double wtable[], int offset)
     {
          final double twopi=2.0D*Math.PI; //6.28318530717959;
          int     idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
          double  dc2, ai1, ai2, ar1, ar2, ds2;
          int     nbd;
          double  dcp, arg, dsp, ar1h, ar2h;
          int iw1 = offset;
          
          arg=twopi / (double)ip;
          dcp=Math.cos(arg);
          dsp=Math.sin(arg);
          nbd=(ido-1)/ 2;
          ipph=(ip+1)/ 2;
          if(ido>=l1)
          {
	      for(k=0; k<l1; k++)
	      {
	          for(i=0; i<ido; i++)
	          {
		      ch[i+k*ido]=cc[i+k*ip*ido];
	          }
	      }
          }
          else
          {
	      for(i=0; i<ido; i++)
	      {
	          for(k=0; k<l1; k++)
	          {
		      ch[i+k*ido]=cc[i+k*ip*ido];
	          }
	      }
          }
          for(j=1; j<ipph; j++)
          {
	      jc=ip-j;
	      j2=2*j;
	      for(k=0; k<l1; k++)
	      {
	          ch[(k+j*l1)*ido]=cc[ido-1+(j2-1+k*ip)*ido]+cc[ido-1+(j2-1+k*ip)*ido];
	          ch[(k+jc*l1)*ido]=cc[(j2+k*ip)*ido]+cc[(j2+k*ip)*ido];
	      }
          }

          if(ido !=1)
          {
	      if(nbd>=l1)
	      {
	          for(j=1; j<ipph; j++)
	          {
		      jc=ip-j;
		      for(k=0; k<l1; k++)
		      {
		          for(i=2; i<ido; i+=2)
		          {
			      ic=ido-i;
			      ch[i-1+(k+j*l1)*ido]=cc[i-1+(2*j+k*ip)*ido]+cc[ic-1+(2*j-1+k*ip)*ido];
			      ch[i-1+(k+jc*l1)*ido]=cc[i-1+(2*j+k*ip)*ido]-cc[ic-1+(2*j-1+k*ip)*ido];
			      ch[i+(k+j*l1)*ido]=cc[i+(2*j+k*ip)*ido]-cc[ic+(2*j-1+k*ip)*ido];
			      ch[i+(k+jc*l1)*ido]=cc[i+(2*j+k*ip)*ido]+cc[ic+(2*j-1+k*ip)*ido];
		          }
		      }
	          }
	      }
	      else
	      {
	          for(j=1; j<ipph; j++)
	          {
		      jc=ip-j;
		      for(i=2; i<ido; i+=2)
		      {
		          ic=ido-i;
		          for(k=0; k<l1; k++)
		          {
			      ch[i-1+(k+j*l1)*ido]=cc[i-1+(2*j+k*ip)*ido]+cc[ic-1+(2*j-1+k*ip)*ido];
			      ch[i-1+(k+jc*l1)*ido]=cc[i-1+(2*j+k*ip)*ido]-cc[ic-1+(2*j-1+k*ip)*ido];
			      ch[i+(k+j*l1)*ido]=cc[i+(2*j+k*ip)*ido]-cc[ic+(2*j-1+k*ip)*ido];
			      ch[i+(k+jc*l1)*ido]=cc[i+(2*j+k*ip)*ido]+cc[ic+(2*j-1+k*ip)*ido];
		          }
		      }
	          }
	      }
          }

          ar1=1;
          ai1=0;
          for(l=1; l<ipph; l++)
          {
	      lc=ip-l;
	      ar1h=dcp*ar1-dsp*ai1;
	      ai1=dcp*ai1+dsp*ar1;
	      ar1=ar1h;
	      for(ik=0; ik<idl1; ik++)
	      {
	         c2[ik+l*idl1]=ch2[ik]+ar1*ch2[ik+idl1];
	         c2[ik+lc*idl1]=ai1*ch2[ik+(ip-1)*idl1];
	      }
	      dc2=ar1;
	      ds2=ai1;
	      ar2=ar1;
	      ai2=ai1;
	      for(j=2; j<ipph; j++)
	      {
	         jc=ip-j;
	         ar2h=dc2*ar2-ds2*ai2;
	         ai2=dc2*ai2+ds2*ar2;
	         ar2=ar2h;
	         for(ik=0; ik<idl1; ik++)
	         {
		     c2[ik+l*idl1]+=ar2*ch2[ik+j*idl1];
		     c2[ik+lc*idl1]+=ai2*ch2[ik+jc*idl1];
	         }
	      }
          }
          for(j=1; j<ipph; j++)
          {
	      for(ik=0; ik<idl1; ik++)
	      {
	         ch2[ik]+=ch2[ik+j*idl1];
	      }
          }
          for(j=1; j<ipph; j++)
          {
	      jc=ip-j;
	      for(k=0; k<l1; k++)
	      {
	          ch[(k+j*l1)*ido]=c1[(k+j*l1)*ido]-c1[(k+jc*l1)*ido];
	          ch[(k+jc*l1)*ido]=c1[(k+j*l1)*ido]+c1[(k+jc*l1)*ido];
	      }
          }

          if(ido==1) return;
          if(nbd>=l1)
          {
	      for(j=1; j<ipph; j++)
	      {
	          jc=ip-j;
	          for(k=0; k<l1; k++)
	          {
		      for(i=2; i<ido; i+=2)
		      {
		          ch[i-1+(k+j*l1)*ido]=c1[i-1+(k+j*l1)*ido]-c1[i+(k+jc*l1)*ido];
		          ch[i-1+(k+jc*l1)*ido]=c1[i-1+(k+j*l1)*ido]+c1[i+(k+jc*l1)*ido];
		          ch[i+(k+j*l1)*ido]=c1[i+(k+j*l1)*ido]+c1[i-1+(k+jc*l1)*ido];
		          ch[i+(k+jc*l1)*ido]=c1[i+(k+j*l1)*ido]-c1[i-1+(k+jc*l1)*ido];
		      }
	          }
	      }
          }
          else
          {
	      for(j=1; j<ipph; j++)
	      {
	          jc=ip-j;
	          for(i=2; i<ido; i+=2)
	          {
		      for(k=0; k<l1; k++)
		      {
		          ch[i-1+(k+j*l1)*ido]=c1[i-1+(k+j*l1)*ido]-c1[i+(k+jc*l1)*ido];
		          ch[i-1+(k+jc*l1)*ido]=c1[i-1+(k+j*l1)*ido]+c1[i+(k+jc*l1)*ido];
		          ch[i+(k+j*l1)*ido]=c1[i+(k+j*l1)*ido]+c1[i-1+(k+jc*l1)*ido];
		          ch[i+(k+jc*l1)*ido]=c1[i+(k+j*l1)*ido]-c1[i-1+(k+jc*l1)*ido];
		      }
	          }
	      }
          }
          for(ik=0; ik<idl1; ik++) c2[ik]=ch2[ik];
          for(j=1; j<ip; j++)
	     for(k=0; k<l1; k++)
	        c1[(k+j*l1)*ido]=ch[(k+j*l1)*ido];
          if(nbd<=l1)
          {
	     is=-ido;
	     for(j=1; j<ip; j++)
	     {
	         is+=ido;
	         idij=is-1;
	         for(i=2; i<ido; i+=2)
	         {
		     idij+=2;
		     for(k=0; k<l1; k++)
		     {
		         c1[i-1+(k+j*l1)*ido] = wtable[idij-1+iw1]*ch[i-1+(k+j*l1)*ido]
                                               -wtable[idij+iw1]*ch[i+(k+j*l1)*ido];
		         c1[i+(k+j*l1)*ido] = wtable[idij-1+iw1]*ch[i+(k+j*l1)*ido]
                                             +wtable[idij+iw1]*ch[i-1+(k+j*l1)*ido];
		     }
	         }
	     }
          }
          else
          {
	     is=-ido;
	     for(j=1; j<ip; j++)
	     {
	         is+=ido;
	         for(k=0; k<l1; k++)
	         {
		     idij=is-1;
		     for(i=2; i<ido; i+=2)
		     {
		         idij+=2;
		         c1[i-1+(k+j*l1)*ido] = wtable[idij-1+iw1]*ch[i-1+(k+j*l1)*ido]
                                               -wtable[idij+iw1]*ch[i+(k+j*l1)*ido];
		         c1[i+(k+j*l1)*ido] = wtable[idij-1+iw1]*ch[i+(k+j*l1)*ido]
                                             +wtable[idij+iw1]*ch[i-1+(k+j*l1)*ido];
		     }
	         }
	     }
          }
     } 

/*---------------------------------------------------------
   rfftf1: further processing of Real forward FFT
  --------------------------------------------------------*/
     void rfftf1(int n, double c[], final double wtable[], int offset)
     {
          int     i;
          int     k1, l1, l2, na, kh, nf, ip, iw, ido, idl1;

          double[] ch = new double[n];
          System.arraycopy(wtable, offset, ch, 0, n);

          nf=(int)wtable[1+2*n+offset];
          na=1;
          l2=n;
          iw=n-1+n+offset;
          for(k1=1; k1<=nf;++k1)
          {
	      kh=nf-k1;
	      ip=(int)wtable[kh+2+2*n+offset];
	      l1=l2 / ip;
	      ido=n / l2;
	      idl1=ido*l1;
	      iw-=(ip-1)*ido;
	      na=1-na;
	      if(ip==4)
	      {
	          if(na==0)
                  {
                      radf4(ido, l1, c, ch, wtable, iw);
                  }
	          else
                  {
                      radf4(ido, l1, ch, c, wtable, iw); 
                  }
	      }
	      else if(ip==2)
	      {
	          if(na==0)
                  {
                      radf2(ido, l1, c, ch, wtable, iw);
                  }
	          else
                  {
                      radf2(ido, l1, ch, c, wtable, iw);
                  }
	      }
	      else if(ip==3)
	      {
	          if(na==0)
                  {
                        radf3(ido, l1, c, ch, wtable, iw);
                  }
	          else
                  {
                        radf3(ido, l1, ch, c, wtable, iw);
                  }
	      }
	      else if(ip==5)
	      {
	          if(na==0)
                  {
                       radf5(ido, l1, c, ch, wtable, iw);
                  }
	          else
                  {
                       radf5(ido, l1, ch, c, wtable, iw);
                  }
	      }
	      else
	      {
                  if(ido==1) na=1-na;
	          if(na==0)
	          {
                      radfg(ido, ip, l1, idl1, c, c, c, ch, ch, wtable, iw);
		      na=1;
	          }
	          else
	          {
                      radfg(ido, ip, l1, idl1, ch, ch, ch, c, c, wtable, iw);
		      na=0;
	          }
	      }
	      l2=l1;
          }
          if(na==1) return;
          for(i=0; i<n; i++) c[i]=ch[i];
     }

/*---------------------------------------------------------
   rfftb1: further processing of Real backward FFT
  --------------------------------------------------------*/
     void rfftb1(int n, double c[], final double wtable[], int offset)
     {
          int     i;
          int     k1, l1, l2, na, nf, ip, iw, ido, idl1;

          double[] ch = new double[n];
          System.arraycopy(wtable, offset, ch, 0, n);

          nf=(int)wtable[1+2*n+offset];
          na=0;
          l1=1;
          iw=n+offset;
          for(k1=1; k1<=nf; k1++)
          {
	      ip=(int)wtable[k1+1+2*n+offset];
	      l2=ip*l1;
	      ido=n / l2;
	      idl1=ido*l1;
	      if(ip==4)
	      {
	          if(na==0) 
                  {
                          radb4(ido, l1, c, ch, wtable, iw);
                  }
	          else
                  {
                          radb4(ido, l1, ch, c, wtable, iw);
                  }
	          na=1-na;
	      }
	      else if(ip==2)
	      {
	          if(na==0)
                  {
                        radb2(ido, l1, c, ch, wtable, iw);
                  }
	          else
                  {
                        radb2(ido, l1, ch, c, wtable, iw);
                  }
	          na=1-na;
	      }
	      else if(ip==3)
	      {
	          if(na==0)
                  {
                          radb3(ido, l1, c, ch, wtable, iw);
                  }
	          else
                  {
                          radb3(ido, l1, ch, c, wtable, iw);
                  }
	          na=1-na;
	      }
	      else if(ip==5)
	      {
	          if(na==0)
                  {
                          radb5(ido, l1, c, ch, wtable, iw);
                  }
	          else
                  {
                          radb5(ido, l1, ch, c, wtable, iw);
                  }
	          na=1-na;
	      }
	      else
	      {
	          if(na==0)
                  {
                          radbg(ido, ip, l1, idl1, c, c, c, ch, ch, wtable, iw);
                  }
	          else
                  {
                          radbg(ido, ip, l1, idl1, ch, ch, ch, c, c, wtable, iw);
                  }
	          if(ido==1) na=1-na;
	      }
	      l1=l2;
	      iw+=(ip-1)*ido;
         }
         if(na==0) return;
         for(i=0; i<n; i++) c[i]=ch[i];
     }

/*---------------------------------------------------------
   rfftf: Real forward FFT
  --------------------------------------------------------*/
     void rfftf(int n, double r[], double wtable[])
     {
         if(n==1) return;
         rfftf1(n, r, wtable, 0);
     } 	/*rfftf*/

/*---------------------------------------------------------
   rfftf: Real backward FFT
  --------------------------------------------------------*/
     void rfftb(int n, double r[], double wtable[])
     {
         if(n==1) return;
         rfftb1(n, r, wtable, 0);
     } /*rfftb*/

/*---------------------------------------------------------
   rffti1: further initialization of Real FFT
  --------------------------------------------------------*/
     void rffti1(int n, double wtable[], int offset)
     {

         final int[] ntryh= new int[] {4, 2, 3, 5};
         final double twopi=2.0D*Math.PI;
         double  argh;
         int     ntry=0, i, j;
         double  argld;
         int     k1, l1, l2, ib;
         double  fi;
         int     ld, ii, nf, ip, nl, is, nq, nr;
         double  arg;
         int     ido, ipm;
         int     nfm1;

         nl=n;
         nf=0;
         j=0;

       factorize_loop:
         while(true)
         {
             ++j;
             if(j<=4)
	         ntry=ntryh[j-1];
             else
	         ntry+=2;
             do
             {
	         nq=nl / ntry;
	         nr=nl-ntry*nq;
	         if(nr !=0) continue factorize_loop;
	         ++nf;
                 wtable[nf+1+2*n+offset]=ntry;

	         nl=nq;
	         if(ntry==2 && nf !=1)
	         {
	             for(i=2; i<=nf; i++)
	             {
		         ib=nf-i+2;
                         wtable[ib+1+2*n+offset]=wtable[ib+2*n+offset];
	             }
                     wtable[2+2*n+offset]=2;
	         }
             }while(nl !=1);
             break factorize_loop;
         }
         wtable[0+2*n+offset] = n;
         wtable[1+2*n+offset] = nf;
         argh=twopi /(double)(n);
         is=0;
         nfm1=nf-1;
         l1=1;
         if(nfm1==0) return;
         for(k1=1; k1<=nfm1; k1++)
         {
             ip=(int)wtable[k1+1+2*n+offset];
	     ld=0;
	     l2=l1*ip;
	     ido=n / l2;
	     ipm=ip-1;
	     for(j=1; j<=ipm;++j)
	     {
	         ld+=l1;
	         i=is;
	         argld=(double)ld*argh;

	         fi=0;
	         for(ii=3; ii<=ido; ii+=2)
	         {
		     i+=2;
		     fi+=1;
		     arg=fi*argld;
                     wtable[i-2+n+offset] = Math.cos(arg);
                     wtable[i-1+n+offset] = Math.sin(arg);
	         }
	         is+=ido;
	     }
	     l1=l2;
         }
     } /*rffti1*/

/*---------------------------------------------------------
   rffti: Initialization of Real FFT
  --------------------------------------------------------*/
     void rffti(int n, double wtable[])  /* length of wtable = 2*n + 15 */
     {
         if(n==1) return;
         rffti1(n, wtable, 0);
     } /*rffti*/

}
